#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix crossDeltaMLOG(NumericMatrix MM,NumericMatrix HH,NumericVector A) {
 NumericMatrix CD(MM.nrow(),MM.ncol());
  for(int m=0;m<MM.nrow();m++) {
    for(int j=0;j<MM.ncol();j++) {
      if(MM(m,j)==1) {
        int N = HH.nrow();
        double total = 0;
        for(int h = 0; h < N; ++h) {
          total += A(h)*HH(h,m)*HH(h,j);
        }
        CD(m,j) = total/N;
      }
    }  
  }
  return CD;
}

// [[Rcpp::export]]
NumericMatrix crossElasMLOG(NumericMatrix MM,NumericMatrix HH,NumericVector A,NumericVector S, NumericVector P) {
  NumericMatrix CD(MM.nrow(),MM.ncol());
  for(int m=0;m<MM.nrow();m++) {
    for(int j=0;j<MM.ncol();j++) {
        int N = HH.nrow();
        double total = 0;
        for(int h = 0; h < N; ++h) {
          total += A(h)*HH(h,m)*HH(h,j)*P(m)/S(j);
        }
        CD(m,j) = total/N;
    }  
  }
  return CD;
}




// [[Rcpp::export]]
NumericMatrix crossDeltaMCES(NumericMatrix MM,NumericMatrix HH,NumericVector A, NumericVector S, NumericVector y) {
  NumericMatrix CD(MM.nrow(),MM.ncol());
  for(int m=0;m<MM.nrow();m++) {
    for(int j=0;j<MM.ncol();j++) {
      if(MM(m,j)==1) {
        int N = HH.nrow();
        double total = 0;
        for(int h = 0; h < N; ++h) {
          total += A(h)*HH(h,m)*HH(h,j)*y(h)/S(j);
        }
        CD(m,j) = total;
      }
    }  
  }
  return CD;
}


// [[Rcpp::export]]
NumericMatrix crossElasMCES(NumericMatrix MM,NumericMatrix HH,NumericVector A, NumericVector S, NumericVector y) {
  NumericMatrix CD(MM.nrow(),MM.ncol());
  for(int m=0;m<MM.nrow();m++) {
    for(int j=0;j<MM.ncol();j++) {
        int N = HH.nrow();
        double total = 0;
        for(int h = 0; h < N; ++h) {
          total += A(h)*HH(h,m)*HH(h,j)*y(h)/S(j);
        }
        CD(m,j) = total;
    }  
  }
  return CD;
}

